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ABSTRACT 

This  study  explores  the  assimilation  of  Doppler  radar  radial  velocity  observations  for  cloud-resolving 
hurricane  analysis,  initialization,  and  prediction  with  an  ensemble  Kalman  filter  (EnKF).  The  case  studied 
is  Hurricane  Humberto  (2007),  the  first  landfalling  hurricane  in  the  United  States  since  the  end  of  the  2005  hur¬ 
ricane  season  and  the  most  rapidly  intensifying  near-landfall  storm  in  U.S.  history.  The  storm  caused  extensive 
damage  along  the  southeast  Texas  coast  but  was  poorly  predicted  by  operational  models  and  forecasters.  It  is 
found  that  the  EnKF  analysis,  after  assimilating  radial  velocity  observations  from  three  Weather  Surveillance 
Radars-1988  Doppler  (WSR-88Ds)  along  the  Gulf  coast,  closely  represents  the  best-track  position  and  intensity 
of  Humberto.  Deterministic  forecasts  initialized  from  the  EnKF  analysis,  despite  displaying  considerable 
variability  with  different  lead  times,  are  also  capable  of  predicting  the  rapid  formation  and  intensification  of 
the  hurricane.  These  forecasts  are  also  superior  to  simulations  without  radar  data  assimilation  or  with  a  three- 
dimensional  variational  scheme  assimilating  the  same  radar  observations.  Moreover,  nearly  all  members  from 
the  ensemble  forecasts  initialized  with  EnKF  analysis  perturbations  predict  rapid  formation  and  intensification 
of  the  storm.  However,  the  large  ensemble  spread  of  peak  intensity,  which  ranges  from  a  tropical  storm  to  a 
category  2  hurricane,  echoes  limited  predictability  in  deterministic  forecasts  of  the  storm  and  the  potential  of 
using  ensembles  for  probabilistic  forecasts  of  hurricanes. 


1.  Introduction 

Landfalling  hurricanes  are  among  the  deadliest  and 
costliest  natural  hazards.  Over  the  past  decade,  signifi- 
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cant  progress  has  been  made  in  short-range  track  fore¬ 
casts  of  tropical  cyclones.  The  current-day  average  48-h 
forecast  position  is  as  accurate  as  a  24-h  track  forecast 
was  10  yr  ago  (Franklin  2004).  However,  there  is  virtu¬ 
ally  no  improvement  in  our  ability  to  predict  hurricane 
intensity  in  terms  of  minimum  sea  level  pressure,  max¬ 
imum  wind  speed,  or  amount  of  precipitation  (Houze 
et  al.  2007).  We  thus  have  very  limited  skill  in  pre¬ 
dicting  tropical  cyclone  formation,  rapid  intensification. 
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fluctuation,  or  decay  (Elsberry  et  al.  2007).  High-resolution 
cloud-resolving  mesoscale  models,  along  with  better 
initialization  of  the  initial  vortex,  may  be  necessary  to 
faithfully  represent  the  internal  dynamics  that  is  crucial 
for  hurricane  intensity  forecasts  (Houze  et  al.  2007; 
Chen  et  al.  2007;  Davis  et  al.  2008). 

Despite  improvements  in  using  advanced  data  assimi¬ 
lation  methods  with  or  without  initial  vortex  bogussing,  our 
ability  to  initialize  a  tropical  cyclone  with  dynamically 
consistent  structure  and  intensity  remains  limited,  even 
with  the  assimilation  of  radar  observations  (e.g.,  Zou  and 
Xiao  2000;  Pu  and  Braun  2001;  Xiao  et  al.  2007).  Nu¬ 
merical  weather  prediction  models  also  have  known 
difficulties  in  their  “spinup”  of  a  tropical  cyclone  or  hurri¬ 
cane  vortex  with  appropriate  moisture,  diabatic,  and  di¬ 
vergence  structures  at  the  initial  time.  Part  of  the  difficulty 
of  hurricane  initialization  comes  from  the  lack  of  routine 
four-dimensional  observations  with  sufficient  spatial  and 
temporal  resolution  to  represent  the  initial  hurricane 
structure  and  intensity.  Another  part  of  the  difficulty  co¬ 
mes  from  the  deficiency  of  the  current  generation  of  op¬ 
erational  data  assimilation  systems,  which  use  static 
background  error  covariance.  The  mostly  balanced,  iso¬ 
tropic,  flow-independent  background  statistics  derived 
from  long-term  averages  of  past  short-term  forecast  error 
(Parrish  and  Derber  1992)  are  ill-suited  for  the  highly 
flow-dependent  background  error  covariances  associated 
with  tropical  cyclones.  In  addition,  operational  models 
generally  have  insufficient  model  resolution  to  effectively 
incorporate  high-resolution  convective-scale  observations 
(such  as  those  from  radars)  for  cloud-resolving  hurricane 
prediction.  Physical  (diabatic)  initializations  using  rainfall, 
radar,  and/or  satellite  observations  are  a  promising  ap¬ 
proach  (Krishnamurti  et  al.  2001),  though  its  effectiveness 
in  spinning  up  a  full  hurricane  vortex  for  cloud-resolving 
hurricane  prediction  remains  to  be  fully  explored. 

The  ensemble  Kalman  filter  (EnKF)  is  a  state-estimation 
technique  that  uses  short-term  ensemble  forecasts  to 
estimate  flow-dependent  background  error  covariance 
or  other  probabilistic  aspects  of  the  background  fore¬ 
cast.  It  was  first  proposed  by  Evensen  (1994)  and  has 
been  adopted  for  data  assimilation  of  many  disciplines 
in  the  geosciences  and  beyond  (Evensen  2003;  Hamill 
2006).  For  the  past  few  years,  the  feasibility  and  per¬ 
formance  of  the  EnKF  have  been  demonstrated  through 
both  simulated  and  real-data  observations  ranging  from 
convective  scales  using  radar  observations  (e.g.,  Snyder 
and  Zhang  2003;  Zhang  et  al.  2004;  Dowell  et  al.  2004; 
Tong  and  Xue  2005)  to  mesoscale  and  regional  scales 
(e.g.,  Zhang  et  al.  2006a;  Torn  et  al.  2006;  Meng  and 
Zhang  2007;  Fujita  et  al.  2007;  Meng  and  Zhang  2008a,b). 
The  benefits  of  directly  estimating  forecast  covariances 
are  also  likely  to  increase  with  the  next  generation  of 


NWP  models  that  resolve  scales  at  which  physical  bal¬ 
ances  are  not  amenable  to  stationary,  isotropic  covari¬ 
ance  models  currently  used  for  operational  forecasts. 

The  current  study  explores  for  the  first  time  the  use 
of  EnKF  to  directly  assimilate  Doppler  radar  radial 
velocity  observations  for  cloud-resolving  hurricane  anal¬ 
ysis  and  prediction,  both  deterministically  and  prob¬ 
abilistically.  The  case  to  be  examined  is  Hurricane 
Humberto  (2007),  the  first  landfalling  hurricane  in  the 
United  States  since  the  end  of  the  2005  hurricane  season 
and  the  most  rapidly  intensifying  near-landfall  storm  in 
U.S.  history.  Humberto  strengthened  from  a  40  mile  per 
hour  (mph)  depression  at  1200  UTC  12  September  2007 
to  a  92-mph  hurricane  at  0700  UTC  13  September,  a 
52  mph  increase  in  surface  wind  speed  in  only  19  h.  The 
storm  caused  extensive  damage  along  the  southeast  Texas 
coast  and  was  poorly  predicted  by  operational  models  and 
forecasters.  The  real-time  forecast  by  the  operational 
Global  Forecast  System  (GFS)  running  at  the  National 
Centers  for  Environmental  Prediction  (NCEP)  failed  to 
capture  the  intensification  and  genesis  of  the  storm  (Fig.  1). 
The  Weather  Research  and  Forecasting  (WRF)  model 
also  failed  in  postevent,  4.5-km,  cloud-resolving  simula¬ 
tions  that  were  initialized  with  the  GFS  analyses  (as  in  the 
control  ensemble  analysis  and  forecast  detailed  in  subse¬ 
quent  sections)  with  lead  times  every  6  h  from  6  to  48  h. 

At  different  stages  of  formation  and  intensification, 
Humberto  was  within  range  of  coastal  Weather  Sur¬ 
veillance  Radars-1988  Dopplers  (WSR-88Ds)  at  Corpus 
Christ!  (KCRP)  and  Houston-Galveston  (KHGX)  in 
Texas  and  Lake  Charles  (KLCH)  in  Louisiana.  These 
radars  provided  valuable  convective-scale  observations 
of  the  storm  but  were  not  assimilated  by  real-time 
NCEP  operational  models.  This  study  is  among  the  first 
to  apply  the  EnKF  to  assimilate  real-data  radar  velocity 
observations  for  complex  weather  phenomena  exhibit¬ 
ing  different  scales  of  motion  [moving  beyond  a  single 
supercell  storm  examined  in  Dowell  et  al.  (2004)].  The 
following  section  will  introduce  the  forecast  model,  the 
EnKF  technique,  and  the  processing  of  the  observations 
to  be  assimilated.  Section  3  presents  the  use  of  the 
EnKF  for  this  storm  in  terms  of  analysis  quality,  and 
deterministic  and  ensemble  forecasts.  Comparison  to  the 
performance  of  data  assimilation  with  a  3-dimensional 
variational  method  that  assimilates  the  same  radar  ob¬ 
servations  as  in  the  WRF  model  is  given  in  section  4. 
Concluding  remarks  are  given  in  section  5. 

2.  Methodology 

a.  The  forecast  model:  WRF 

The  Advanced  Research  WRF  (ARW)  is  used  in 
this  study.  WRF  is  a  fully  compressible,  nonhydrostatic 
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Fig.  1.  Time  evolution  of  (a)  minSLP  and  (b)  maxWSP  forecasts  by  operational  GFS  forecasts  starting  every  6  h 
from  0000  UTC  12  Sep  to  0000  UTC  12  Sep  2007  and  by  4.5-km  WRF  forecasts  starting  from  the  operational  GFS 
analyses  in  comparison  to  the  NHC  best-track  estimate. 


mesoscale  model  (Skamarock  et  al.  2005).  The  vertical 
coordinate  follows  the  terrain  using  hydrostatic  pres¬ 
sure,  and  the  model  uses  an  Arakawa  C  grid.  Prognostic 
variables  are  the  column  mass  of  dry  air,  velocities  (w,  v,  w), 
potential  temperature,  geopotential,  and  mixing  ratios 
for  water  vapor,  cloud,  rain,  ice,  snow,  and  graupel. 

In  the  control  experiments,  three  model  domains  with 
two-way  nesting  are  used  (Fig.  2).  The  two  coarse  do¬ 
mains  (D1  and  D2)  both  have  160  X  121  grid  points  and 
grid  spacings  of  40.5  and  13.5  km,  respectively.  The  in¬ 
nermost  domain,  D3,  has  253  X  253  grid  points  and  a 
grid  spacing  of  4.5  km.  All  model  domains  have  35 
vertical  layers,  and  the  model  top  is  set  at  10  hPa.  The 
physical  parameterization  schemes  include  the  Grell- 
Devenyi  cumulus  scheme  (Grell  and  Devenyi  2002), 
WRF  single-moment  six-class  microphysics  with  grau¬ 
pel  (Hong  et  al.  2004),  and  the  Yonsei  State  University 
(YSU)  scheme  (Noh  et  al.  2003)  for  planetary  boundary 
layer  processes.  The  NCEP  GFS  operational  analysis  at 
0000  UTC  12  September  and  its  forecast  are  used  to 
create  the  initial  and  boundary  conditions.  Data  as¬ 
similation  is  performed  for  all  domains  but  all  verifica¬ 
tion  is  performed  for  D3. 

b.  The  data  assimilation  method:  EnKF 

The  EnKF  implemented  in  the  WRF  model  is  the 
same  as  that  in  Meng  and  Zhang  (2008a,b)  except  that 
no  multischeme  ensemble  is  used  for  this  study.  This 
version  of  the  filter  was  originally  implemented  in  the 
fifth-generation  Pennsylvania  State  University-National 
Center  for  Atmospheric  Research  Mesoscale  Model 
(MM5),  which  is  documented  in  Zhang  et  al.  (2006a).  It 
uses  the  covariance  relaxation  of  Zhang  et  al.  (2004)  to 
inflate  the  background  error  covariance.  Unlike  the 
standard  inflation  method  (Anderson  2001),  in  which  all 


points  in  the  prior  field  are  inflated,  this  relaxation 
method  only  inflates  the  covariance  at  updated  points 
via  a  weighted  average  between  the  prior  perturbation 
(denoted  by  superscript  /)  and  the  posterior  perturba¬ 
tion  (denoted  by  superscript  a)  as  follows: 

«,J'  =  (l-a)(x“)'+a(x0'.  (1) 

The  weighting  coefficient,  a,  is  set  to  0.5  in  the  observing 
system  simulation  experiment  (OSSE)  studies  of  Zhang 
et  al.  (2004,  2006a)  for  the  perfect-model  experiments 
but  a  range  from  0.7  to  0.8  is  found  to  be  necessary  for 
imperfect-model  experiments  (Meng  and  Zhang  2007) 
or  real-data  applications  (Whitaker  et  al.  2008;  Meng 
and  Zhang  2008a,b;  Torn  and  Hakim  2008)  due  to  un¬ 
avoidable  imperfections  in  the  forecast  model.  A  value 
of  0.8  is  used  for  the  present  real-data  study. 

c.  Ensemble  initial  and  boundary  perturbations 

Although  the  optimum  ensemble  size  for  estimating 
the  forecast  uncertainty  is  still  under  active  research,  30 
members  are  used  herein.  An  ensemble  size  of  20-50  was 
found  to  be  affordable  and  reasonable  based  on  previous 
studies  (e.g.,  Houtekamer  and  Mitchell  2001;  Anderson 
2001;  Snyder  and  Zhang  2003;  Zhang  2005;  Zhang 
et  al.  2006a;  Meng  and  Zhang  2007;  Meng  and  Zhang 
2008a,b).  As  in  Zhang  et  al.  (2006a),  the  initial  ensemble 
is  generated  with  the  WRF’s  three-dimensional  varia¬ 
tional  data  assimilation  (3DVAR)  using  the  cv3  back¬ 
ground  error  covariance  option  (Barker  et  al.  2004).  To 
create  a  largely  balanced  perturbation,  we  first  generate 
a  set  of  random  control  vectors  with  a  normal  distri¬ 
bution  (zero  mean  and  unit  standard  deviation).  Then, 
the  control  increment  vector  is  transformed  back  to 
model  space  via  an  EOF  transform,  a  recursive  filter. 
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Fig.  2.  Configuration  of  WRF  model  domains  1, 2,  and  3  with  horizontal  grid  spacings  of  40.5, 
13.5,  and  4.5  km,  respectively.  Also  depicted  are  the  NHC  best-track  estimates  of  Humberto 
with  intensity  in  gray  scale  and  the  three  WSR-88D  locations. 


and  physical  transformation  via  balance  equations.  The 
perturbed  variables  include  horizontal  wind  compo¬ 
nents,  potential  temperature,  and  mixing  ratio  for  water 
vapor,  and  their  error  statistics  are  defined  by  the  cli¬ 
matological  background  error  covariance.  Other  prog¬ 
nostic  variables  such  as  vertical  velocity  (iv)  and  mixing 
ratios  for  cloud  water  rainwater  snow  (^^),  and 
graupel  {qg)  are  not  perturbed.  The  perturbation  stan¬ 
dard  deviations  thus  generated  are  approximately 
2  m  s~^  for  horizontal  wind  components  {u  and  v),  0.8  K 
for  temperature  (T),  1  hPa  for  pressure  perturbation 
(p'),  and  0.8  g  kg“^  for  the  water  vapor  mixing  ratio  (q). 
The  3DVAR  perturbations  are  added  to  the  GFS 
analysis  to  form  an  initial  ensemble,  which  is  then  in¬ 
tegrated  for  9  h  to  develop  an  approximately  realistic, 
flow-dependent  background  error  covariance  structure 
before  the  first  observation  is  assimilated.  Similar  meth¬ 
ods,  using  3DVAR  to  generate  the  initial  ensemble  for 
the  EnKF,  are  also  employed  in  Houtekamer  et  al. 
(2005)  and  Barker  (2005). 

The  simplest  way  to  perturb  lateral  boundary  condi¬ 
tions  for  a  limited-area  model  is  to  use  a  global  en¬ 
semble  forecast  with  the  correct  size  and  resolution 
[which  are  usually  unavailable;  Chessa  et  al.  (2004)]. 
Torn  et  al.  (2006)  examined  several  alternative  bound¬ 
ary  perturbation  methods  and  concluded  that  the  error 
originating  from  using  different  methods  is  limited  to 
near  the  edges  of  the  domain.  In  this  paper,  the  GFS 


operational  forecasts  are  used  to  create  boundary  con¬ 
ditions  that  are  perturbed  in  the  same  manner  as  with 
the  initial  ensemble. 

d.  Super  observations  and  quality  control 

With  large  volumes  of  radar  observations  that  are 
recorded  at  a  much  higher  resolution  than  the  forecast 
model  grid  spacing  for  the  EnKF  data  assimilation, 
significant  data  thinning  and  quality  control  of  obser¬ 
vations  become  necessary.  The  process  of  combining 
multiple  observations  into  one  high-accuracy  “super” 
observation  (SO)  is  often  referred  to  as  “superobbing.” 
An  SO  for  radar  radial  velocity  is  created  through 
horizontal  averaging  in  polar  space  of  the  raw  polar 
volume  of  data  (Lindskog  et  al.  2000,  2004;  Alpert  and 
Kumar  2007).  To  minimize  horizontal  correlations  of 
the  SOs,  each  pixel  of  the  raw  data  is  allowed  to  influ¬ 
ence  one  SO  only.  To  avoid  the  averaging  of  the  radial 
velocity  (Vr)  in  significantly  different  directions,  the 
averaging  bin  is  confined  within  5  km  in  the  radial  di¬ 
rection  and  5°  in  the  azimuthal  direction.  Our  selection 
of  the  bin  size  is  within  the  range  of  values  used  in  the 
literature  (e.g.,  Lindskog  et  al.  2000;  Alpert  and  Kumar 
2007;  Montmerle  and  Faccani  2009). 

For  quality  control  of  the  observations,  we  first  use 
the  NCAR  radar  editing  software  called  SOLO  [infor¬ 
mation  online  at  http://www.eol.ucar.edu/rdp/solo/solo_ 
home.html]  to  dealias  the  range-folded  data  and  to 
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remove  apparently  erroneous  observations.  We  then 
implement  the  following  additional  quality  control  mea¬ 
sures  in  the  SO  generation  for  this  study:  1)  any  raw 
observations  with  values  smaller  than  2  m  s~^  or  larger 
than  70  m  s“^  or  with  distances  to  the  radar  smaller  than 
4  km  will  be  discounted,  2)  a  raw  Vr  observation  will  be 
discounted  if  its  deviation  from  the  bin  mean  exceeds 
twice  the  standard  deviation  of  all  raw  observations  in 
the  bin,  3)  there  shall  be  at  least  four  valid  Vr  raw  ob¬ 
servations  within  an  averaging  bin,  4)  there  will  be  no 
SO  for  a  bin  whose  standard  deviation  is  twice  the  av¬ 
erage  of  the  standard  deviations  in  all  bins,  and  5)  the 
final  SO  value  of  the  bin  will  be  the  average  of  at  most 
10  raw  observations  that  are  closest  to  the  center  of  the 
bin.  These  quality  control  procedures  are  in  place  to 
minimize  the  impact  of  ground  clutter  and  to  correct  the 
failures  in  the  subjective  dealiasing  step;  similar  proce¬ 
dures  were  also  employed  in  Xu  et  al.  (2003)  and 
Montmerle  and  Faccani  (2009). 

Additional  quality  controls  are  implemented  in  the 
processing  of  EnKF  analyses.  The  observation  errors  of 
all  SOs  are  assumed  to  be  3  m  s“^  in  this  study,  and  an  SO 
will  also  be  discounted  if  the  difference  between  this  SO 
and  the  forecast  prior  is  larger  than  5  times  the  obser¬ 
vation  error.  The  3  m  s~^  observational  error  is  consistent 
with  the  range  used  in  Dowell  et  al.  (2004)  and  Montmerle 
and  Faccani  (2009),  and  is  also  a  conservative  version  of 
the  2.4  m  s“^  estimated  in  Xu  et  al.  (2003). 

e.  Successive  covariance  localization 

A  successive  covariance  localization  (SCL)  technique 
is  designed  to  assimilate  dense  radar  observations  that 
contain  information  about  the  state  of  the  atmosphere 
at  a  wide  range  of  scales.  The  method  is  also  designed  to 
reduce  computational  costs  and  sampling  errors.  This 
technique  uses  the  Gaspari  and  Cohn  (1999)  fifth-order 
correlation  function  for  covariance  localization,  but  a 
different  localization  radius  of  influence  (ROI)  is  used 
for  different  groups  of  observations  by  random  sam¬ 
pling.  SCL  assumes  that  both  large-  and  small-scale 
errors  are  simultaneously  present.  First,  one  tries  to 
remove  dynamically  important  aspects  of  the  large- 
scale  error  by  assimilating  a  relatively  small  subset  of 
observations  with  a  large  ROI.  Next,  the  ROI  is  made 
smaller,  and  higher-density  observations  are  used  to 
constrain  both  smaller-scale  errors  and  what  remains  of 
the  large-scale  error.  The  process  is  repeated  until  all 
scales  resolved  by  the  observational  network  have  been 
adequately  dealt  with.  The  SCL  method  has  some  re¬ 
semblance  to  the  successive  correction  method  used 
in  some  earlier  empirical  objective  analysis  schemes 
(e.g.,  Barnes  1964),  though  in  the  EnKF  the  same  ob¬ 
servation  will  not  be  used  twice.  Sensitivity  experiments 


that  demonstrate  the  benefits  of  using  the  SCL  method 
over  using  single  ROIs  for  all  observations  will  be  pre¬ 
sented  in  section  3e.  The  use  of  SCL  is  partially  moti¬ 
vated  by  the  fact  that  with  serial  observation  processing 
of  the  EnKF,  the  error  correlation  length  scale  de¬ 
creases  as  the  previously  assimilated  observations  better 
define  the  large  scales;  hence,  later  observations  should 
be  assimilated  with  tighter  localization  (Bishop  and 
Hodyss  2007). 

In  this  particular  case,  we  first  use  a  horizontal  ROI  of 
1215  km  for  10%  of  all  SOs  to  capture  the  large-scale 
background  flow  in  all  three  domains  (the  numbers  of 
SOs  at  each  time  from  each  radar  are  shown  in  Fig.  3). 
We  then  use  an  ROI  of  405  km  to  assimilate  another 
20%  of  the  total  SOs  to  represent  the  mesoscale  flow 
(i.e.,  tropical  cyclone  scale)  for  the  13.5-  and  4.5-km 
domains  (i.e.,  D2  and  D3).  Last,  we  use  an  ROI  of  135  km 
to  assimilate  another  60%  of  the  total  SOs  to  capture 
even  smaller-scale  phenomena  that  include  mesoscale 
vortices  just  in  D3.  In  other  words,  we  assimilate  10%, 
30%,  and  90%  of  the  total  SOs  for  DI,  D2,  and  D3, 
respectively.  While  the  partitioning  of  the  observa¬ 
tions  into  different  bins  remains  subjective,  we  binned 
more  observations  to  smaller  scales  since  arguably 
there  are  larger  degrees  of  freedom  and  lesser  degrees 
of  balance  at  smaller  scales.  The  remaining  10%  of  the 
SOs  that  are  not  assimilated  by  any  domain  will  be  used 
for  verifications. 

The  vertical  ROI  is  34  based  on  the  number  of  ver¬ 
tical  levels  (as  in  Zhang  et  al.  2006a)  for  all  three  do¬ 
mains.  The  Gaspari  and  Cohn  (1999)  fifth-order  corre¬ 
lation  function  is  also  used  for  the  vertical  covariance 
localization. 

3.  EnKF  performance 

a.  EnKF  analyses 

The  control  EnKF  experiment  begins  to  assimilate 
the  SOs  from  the  KCRP  and  KHGX  WSR-88Ds  at  0900 
UTC  12  September  2007,  6  h  before  the  tropical  de¬ 
pression  status  was  declared  by  National  Hurricane 
Center  (NHC).  The  EnKF  continues  to  assimilate  SOs 
from  these  two  radars  every  hour  until  1800  UTC  12 
September.  After  this  time,  the  KCRP  radar  is  too  far 
from  the  storm  to  have  any  significant  impact  while 
KLCH  radar  begins  to  offer  significant  coverage  of  the 
storm.  The  KLCH  radar  observations  are  thus  assimi¬ 
lated  beginning  at  1900  UTC  12  September.  The  EnKF 
assimilation  continues  every  hour  until  1200  UTC  13 
September,  a  few  hours  after  the  storm  in  the  NHC  best 
track  reaches  its  peak  intensity  and  starts  weakening 
over  land.  The  numbers  of  SOs  from  the  three  radars  by 
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Fig.  3.  (a)  The  number  of  SOs  from  each  radar  at  different  times  by  the  control  EnKF 
experiment  and  the  exemplar  distributions  of  SOs  at  (b)  0900  and  (c)  1900  UTC  12  Sep  2007. 


the  EnKF  at  different  times  and  example  distributions 
of  the  SOs  at  0900  and  1900  UTC  12  September  are 
shown  in  Fig.  3. 

Figure  4  shows  the  time  evolution  of  the  minimum  sea 
level  pressure  (minSLP)  and  maximum  surface  wind 
speed  (maxWSP)  estimated  from  the  posterior  EnKF 
mean  analysis  field  (gray;  hereafter  “the  EnKF  anal¬ 
ysis”)  as  well  as  the  average  (red)  of  the  maximum- 
minimum  values  estimated  from  each  ensemble  mem¬ 
ber’s  posterior  (green)  in  comparison  to  the  NHC 
best-track  analyses  (black).  Minimum  SEP  in  the  EnKF 
analysis  (gray)  agrees  well  with  the  average  of  the 
members’  minimum  values  (red)  due  to  collocation  of 
the  centers  of  most  members.  The  analysis  mean  and 
some  ensemble  members  also  compare  favorably  with 
the  NHC  best-track  analysis.  However,  due  to  the 
strong  spatial  and  temporal  variability  of  the  maximum 
surface  wind  speed,  maxWSP  in  the  EnKF  analysis 
(gray)  is  significantly  smaller  than  that  of  the  average 
(red)  of  the  members’  maxima.  Yet,  the  average  of  each 
member’s  maxWSP  (red)  matches  well  with  the  NHC 
best-track  estimate  (black)  but  larger  errors  are  seen 
after  landfall.  We  therefore  believe  it  is  more  appro¬ 


priate  to  use  the  averages  of  maxWSP  and  minSLP  from 
each  member  for  verifying  ensemble  forecasts  in  terms 
of  extreme  values. 

There  is  large  ensemble  spread  of  both  minSLP  and 
maxWSP  among  the  analysis  members  (Fig.  4).  Stan¬ 
dard  deviation  of  minSLP  (maxWSP)  increases  from 
1  to  2  hPa  (1-3  m  s“^)  during  the  first  few  assimilation 
cycles  to  10.2  hPa  (8.3  m  s~^)  near  the  peak  intensity 
time  at  0800  UTC  13  September  but  drops  quickly  to 
5.4  hPa  (4.0  m  s~^)  at  0900  UTC  13  September  after  the 
storm  in  most  members  makes  its  landfall  (not  shown). 
Likewise,  the  minSLP  (maxWSP)  at  the  peak  intensity 
time  varies  from  992  (27  m  s~^)  to  960  hPa  (59  m  s~^).  In 
terms  of  the  corresponding  intensity  categories,  this 
represents  a  range  from  a  strong  tropical  storm  to  a 
major  hurricane.  Large  disparities  between  ensemble 
members  demonstrate  significant  EnKF  analysis  un¬ 
certainties  during  and  after  the  rapid  intensification  of 
Humberto. 

Despite  the  large  spread,  data  assimilation  with  the 
EnKF  is  clearly  beneficial.  All  analysis  ensemble  mem¬ 
bers  capture  the  storm’s  formation  and  intensifica¬ 
tion  when  EnKF  assimilates  Vr  observations.  This  is  in 
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Fig.  4.  Time  evolutions  of  (a)  minSLP  and  (b)  maxWSP  estimated  from  the  EnKF  analysis  and  NoDA  ensemble 
forecast.  Thin  green  (cyan)  lines  represent  the  maxWSP  and  minSLP  estimated  from  each  ensemble  member  of  the 
EnKF  analysis  (NoDA  forecast)  with  thick  red  (blue)  lines  representing  the  average  of  maxWSP  and  minSLP  in  each 
member.  The  gray  line  indicates  the  maxWSP-minSLP  estimated  from  the  EnKF  analysis  mean  while  the  black  curve 
is  the  NHC  best-track  estimate. 


Strong  contrast  to  pure  ensemble  forecasts  started  with 
the  same  prior  perturbations  but  without  EnKF  assim¬ 
ilation  of  Vr  (hereafter  refer  to  as  “NoDA”).  The 
NoDA  ensemble  neither  captures  the  mean  develop¬ 
ment  nor  the  realistic  uncertainties  associated  with  the 
mean  forecast  (blue  curves  in  Fig.  4). 

Figure  5  compares  the  observed  0.5°  base  scan  of  the 
radial  velocity  from  KHGX  (only  1/10  of  the  observa¬ 
tions  plotted)  with  corresponding  simulated  values  from 
both  the  EnKF  analysis  and  the  NoDA  ensemble  fore¬ 
cast  mean  valid  at  0900  and  1800  UTC  12  September, 
and  0300  UTC  13  September.  The  benefit  of  assimilat¬ 
ing  Vr  observations  with  EnKF  is  evident  even  after 
the  first  volume  of  observations  is  assimilated  at  0900 
UTC  12  September.  The  analyses  at  this  time  capture 
the  coastal  mesoscale  circulation  much  better  than  the 
NoDA  ensemble  mean.  Since  NoDA  at  0900  UTC  12 
September  is  simply  the  EnKF  prior  estimate  with  no 
Vr  observations  assimilated,  a  comparison  of  Figs.  5b 
and  5c  (verifying  against  Fig.  5a)  shows  the  immediate 
benefit  of  the  EnKF  at  the  initial  assimilation  time. 
Subsequently,  the  EnKF  well  analyzes  the  cyclone  vortex 
structure,  intensity,  and  evolution  (Fig.  5). 

The  quality  of  the  EnKF  analysis  is  also  indepen¬ 
dently  verified  against  the  additional  10%  of  the  total 
radial  velocity  SOs  from  both  the  KHGX  and  KLCH 
radars  that  have  never  been  assimilated  (refer  to  sec¬ 
tion  2e).  Figure  6  shows  the  root-mean-square  error 
(RMSE)  of  the  EnKF  analysis  in  comparison  to  the 
NoDA  experiment  plotted  every  3  h  (RMSE  is  the  dif¬ 
ference  between  the  simulated  radial  velocity  and 
the  10%  SOs  verified  at  the  SO’s  locations).  Except 
for  at  the  first  analysis  time  for  the  KLCH  radar  and 


consistent  with  the  other  measures  of  analysis  qual¬ 
ity  discussed  above,  the  EnKF  analysis  error  main¬ 
tained  a  remarkably  smaller  amplitude  verifying  against 
these  independent  radial  velocity  SOs.  The  mean  anal¬ 
ysis  errors  verified  against  both  radars  are  similar  to 
the  assumed  observational  error  of  3  m  s“^.  Consistent 
with  the  increasing  analysis  ensemble  spread  in  the 
maxWSP  in  Fig.  3c,  there  is  a  significant  increase  in 
the  analysis  error  during  the  peak  intensity  times,  es¬ 
pecially  when  verifying  against  the  KHGX  radar  SOs 
(Fig.  6a). 

Figure  7  shows  a  comparison  of  the  radar  reflectivities 
from  the  observed  composite,  the  EnKF  analysis,  and 
the  NoDA  ensemble  forecast  mean  valid  at  1200  UTC 
12  September,  and  0000  and  1200  UTC  13  September, 
respectively.  The  impact  of  Vr  observation  assimilation 
on  the  unobserved  reflectivity  variable  is  evident  in  the 
progressively  better  posterior  estimates  (analysis  means) 
of  reflectivity.  At  1200  UTC  12  September,  while  the 
NoDA  ensemble  forecast  mean  simulates  a  broader 
area  of  light  precipitation,  the  EnKF  analysis  begins  to 
localize  the  precipitation  into  two  primary  bands.  The 
first  band  is  located  along  the  Gulf  coast  of  Texas  and 
Louisiana  to  the  east  and  northeast  of  Houston,  and  the 
other  is  farther  west  but  far  south  of  the  display  domain. 
This  compares  much  more  favorably  to  the  observations 
(though  the  convection  in  the  EnKF  analysis  is  still 
weaker  and  broader  partly  due  to  the  ensemble  aver¬ 
aging  effects  discussed  above). 

At  0000  UTC  13  September,  during  the  storm’s  rapid 
intensification,  the  EnKF  analysis  captures  an  im¬ 
pressive  developing  tropical  cyclone  south  of  KHGX 
with  multiple  spiral  rainbands  to  the  north  and  east 
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Fig.  5.  The  raw  radial  velocity  observations  from  the  (left)  0.5°  base  scan  of  the  KHGX  radar  (OBS),  (middle)  corresponding  EnKF 
analysis,  and  (right)  NoDA  ensemble  forecast  mean  valid  at  0900  and  1800  UTC  12  Sep,  and  0300  UTC  13  Sep  2007,  respectively. 


quadrants  fueled  by  warm  moist  air  from  the  south. 
Except  for  a  spurious  onshore  mesoscale  rainband  right 
across  the  Texas  and  Louisiana  border,  the  position, 
intensity,  and  structure  of  the  rainbands,  including  the 
developing  eyewall  in  the  EnKF  analysis  (Fig.  7e), 
compare  remarkably  well  with  the  observed  reflectivity 
(Fig.  7d).  At  1200  UTC  13  September,  the  final  analysis 
time  after  the  storm  begins  its  rapid  weakening  over 
land,  the  EnKF  analysis  correctly  places  the  center  of 
the  storm  over  the  Texas-Louisiana  border.  It  also  cor¬ 


rectly  analyzes  the  broad  rainbands  to  the  east  of  the 
storm  (Figs.  7g  and  7h).  The  NoDA  ensemble  forecast 
mean,  on  the  other  hand,  does  not  simulate  any  tropical 
development  (Figs.  7c,  7f,  and  7i). 

b.  Deterministic  forecasts  from  the  EnKF  analysis 

Next,  we  examine  the  value  that  using  EnKF  to  assim¬ 
ilate  Vr  adds  to  forecasts.  Single,  deterministic  forecasts 
from  the  control  EnKF  mean  analyses  are  performed 
every  3  h  from  1200  UTC  12  September  to  1200  UTC  13 
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Fig.  6.  The  RMSE  of  the  radial  velocity  of  the  EnKF  analysis  vs  the  NoDA  experiment  verified  against  the  10% 
independent  radial  velocity  SOs  from  (a)  KHGX  and  (b)  KLCH. 


September.  All  these  simulations  are  integrated  until 
1200  UTC  14  September. 

Figure  8  shows  the  simulated  cyclone  position, 
minSLP,  and  maxWSP  in  deterministic  WRF  forecasts 
initialized  with  EnKF  analyses  at  different  times.  De¬ 
spite  a  slight  delay  in  peak  intensity  compared  to  the 
best-track  estimates,  all  forecasts  initialized  with  EnKF 
analyses  simulate  significant  tropical  cyclone  formation 
and  intensification  if  WSR-88D  radial  velocity  obser¬ 
vations  are  assimilated  for  9  h  or  longer.  On  average, 
continuous  assimilation  through  time  and  assimilating 
more  observations  produces  both  better  analyses  of  the 
initial  storm  and  better  deterministic  track  and  intensity 
forecasts,  which  is  especially  evident  for  forecasts  ini¬ 
tialized  before  1800  UTC  12  September.  These  WRF 
forecasts  from  the  EnKF  analyses  are  in  strong  contrast 
to  the  nearly  complete  forecast  failure  of  the  WRF 
forecasts  (with  the  same  model  configuration)  cold  star¬ 
ted  from  the  GFS  analyses  at  all  lead  times  (Fig.  1).  This 
signifies  the  importance  and  potential  of  assimilating 
radar  observations  in  improving  cloud-resolving  tropical 
cyclone  initialization  and  prediction  at  all  lead  times. 

Among  the  simulations  initialized  from  the  mean 
EnKF  analyses,  forecasts  from  1800  and  2100  UTC  12 
September  are  the  most  remarkable.  The  peak  intensity 
of  both  simulations  (which  have  significant  lead  times)  is 
within  2  hPa  of  the  observed  (best  track)  peak  intensity. 
Also,  the  maximum  surface  winds  in  both  runs  reach  or 
nearly  reach  category  1  hurricane  intensity,  which  is  less 
than  5  m  s“^  different  from  the  best-track  estimates. 
Both  forecasts  are  considered  to  be  quite  successful 
given  the  uncertainties  in  the  best-track  estimate  and 
since  a  relatively  small  number  of  model  outputs  (every 
3  h)  is  used  for  determining  the  simulated  intensity.  The 
slight  delay  in  the  peak  intensity  for  all  forecasts  reflects 
a  slight  lag  in  the  simulated  landfall  time  compared  to 
the  best-track  observations.  Broadly  speaking,  the  ini¬ 


tial  positions  and  subsequent  track  errors  in  the  EnKF 
analyses  become  smaller  as  the  radar  data  assimilation 
cycle  proceeds. 

Figure  9  compares  the  observed  composite  radar  re¬ 
flectivity  with  the  corresponding  simulated  reflectivity 
in  deterministic  forecasts  valid  at  2100  UTC  12  Sep¬ 
tember  and  0300  and  0900  UTC  13  September.  Exper¬ 
iment  NoDA  is  initialized  at  0900  UTC  12  September 
from  the  ensemble-mean  forecast,  and  experiment 
EnKF  is  initialized  at  1800  UTC  12  September  after  9  h 
of  EnKF  analysis.  Consistent  with  the  track  and  inten¬ 
sity  forecasts  shown  in  Fig.  8,  the  deterministic  forecast 
from  the  EnKF  analysis  mean  compares  favorably  with 
radar  observations  in  terms  of  the  structure  and  place¬ 
ment  of  rainbands  and  the  formation  of  the  eyewall 
right  before  moving  over  the  coast.  These  phenomena 
are  nonexistent  in  the  forecast  without  EnKF  assimila¬ 
tion  of  Vr  observations  (Figs.  9c,  9f,  and  9i).  The  fore¬ 
cast  without  data  assimilation  also  develops  widespread 
convection  over  the  Gulf  region  across  the  display  do¬ 
main,  but  the  convection  is  highly  disorganized  with 
only  weak  cyclonic  circulation  in  the  center.  While  the 
above  forecast  from  the  EnKF  analysis  is  far  better,  it  is 
also  far  from  perfect.  For  example,  the  outer  spiral 
rainband  to  the  far  east  of  the  storm  is  not  well  captured, 
partly  because  this  forecast  is  initialized  at  1800  UTC  12 
September,  before  the  KLCH  radar  observations  are 
assimilated. 

c.  Deterministic  forecast  with  a  1.5-km 
movable  nested  domain 

Since  the  4.5-km  control  experiments  only  marginally 
resolve  moist  convection,  we  perform  another  high- 
resolution  experiment  that  nests  an  additional  domain 
with  1.5-km  horizontal  grid  spacing  (1.5KM).  Experi¬ 
ment  1.5KM  has  a  fourth  domain  of  253  X  253  grid 
points  centered  on  the  maximum  vorticity  and  moves 
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Fig.  7.  Comparison  of  the  (left)  radar  reflectivity  (dBZ)  from  the  observational  composite  mosaic  (OBS),  (middle)  corresponding  EnKF 
analysis,  and  (right)  NoDA  ensemble  forecast  mean  valid  at  1200  UTC  12  Sep,  and  0000  and  1200  UTC  12  Sep  2007,  respectively. 


with  the  tropical  cyclone  using  a  vortex-tracking  method 
(Chen  et  al.  2007)  and  two-way  nesting.  Experiment 
1.5KM  is  initialized  at  1800  UTC  12  September  with  the 
same  EnKF  analysis  from  the  4.5-km  domain  discussed 
above.  Figure  10  shows  the  1.5KM-simulated  minSLP 
and  maxWSP  in  comparison  to  both  the  best-track  es¬ 
timate  and  the  4.5-km  control  experiment  starting  at  the 
same  time.  Aside  from  having  slightly  slower  move¬ 
ment,  the  simulated  track  in  this  1.5-km  forecast  is 
nearly  identical  to  that  of  the  4.5-km  control  forecast 
initialized  at  the  same  time  (Fig.  10a).  The  peak  inten¬ 
sity  is  5  hPa  weaker  in  terms  of  minSLP  (Fig.  10b)  and  is 
comparable  in  terms  of  maxWSP  (Fig.  10c).  The  dif¬ 
ference  is  well  within  the  uncertainties  of  the  EnKF 
analysis  (Fig.  4),  and  the  ensemble  forecasts  started  at 


the  same  time  with  the  EnKF  analysis  perturbations 
(next  subsection). 

Although  the  storm  in  1.5KM  is  slightly  weaker  in 
terms  of  peak  intensity,  Fig.  11  shows  that  the  higher- 
resolution  simulation  captures  a  much  more  realistic 
detailed  mesoscale  structure  (e.g.,  0000  and  0600  UTC 
13  September)  of  the  cyclone.  This  additional  structure 
compares  more  favorably  to  the  observed  radar  re¬ 
flectivity  at  different  stages  of  the  cyclone  development 
than  that  of  the  4.5-km  control  forecast.  In  particular,  it 
is  only  the  1.5KM  run  that  captures  the  strong  asym¬ 
metry  of  the  observed  reflectivity  at  0600  UTC  13 
September.  Future  studies  will  perform  high-resolution 
forecasts  at  different  times  and/or  utilize  the  EnKF 
analysis  on  the  1.5-km  model  domain. 
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Fig.  8.  The  simulated  (a)  positions,  (b)  minSLP,  and  (c) 
maxWSP  of  Humberto  in  the  deterministic  WRF  forecasts  (color 
curves)  initialized  with  the  EnKF  analyses  every  3  h  from  1200 
UTC  12  Sep  to  1200  UTC  13  Sep  2007  in  comparison  with  the 
NHC  best-track  estimate  (black). 


d.  Ensemble  forecasts  from  EnKF  analysis 

Large  variations  between  deterministic  forecasts  ini¬ 
tialized  with  the  EnKF  mean  analysis  at  different  times 
and  the  difficulty  in  real-time  operational  forecasting 
suggest  that  the  predictability  of  this  storm  is  rather 
limited.  The  EnKF  analysis  helps  to  understand  this 
predictability  problem  because  it  provides  consistent, 
flow-dependent  uncertainties  that  are  used  for  initial¬ 
izing  ensembles  for  probabilistic  prediction. 

A  30-h  ensemble  forecast  is  initiated  with  the  control 
EnKF  analysis  and  perturbations  at  1800  UTC  12  Sep¬ 
tember.  Consistent  with  the  large  variations  between 
deterministic  forecasts  starting  with  EnKFs  at  different 
times  (Fig.  8),  there  is  also  large  spread  among  different 
members  from  the  EnKF-initialized  ensemble.  This  is 
shown  in  Fig.  12,  which  plots  the  evolution  of  the  posi¬ 
tion,  minSLP,  and  maxWSP  as  simulated  by  all  mem¬ 
bers  and  the  deterministic  forecast  initialized  from  the 
EnKF  analysis  (mean).  The  spread  of  minSLP  triples 
over  30  h  from  less  than  1.5  hPa  at  1800  UTC  12  Sep¬ 
tember  (the  initial  time)  to  4.5  hPa  at  0900  UTC  13 
September,  while  the  spread  of  maxWSP  grows  from 
1.6  to  6.3  m  s“^  during  the  same  period  (even  larger 
ensemble  spreads  are  observed  at  0700  and  0800  UTC 
13  September;  not  shown). 

The  large  forecast  uncertainty  and  sensitivity  to  initial 
perturbation  can  also  be  clearly  seen  in  Fig.  13,  which 
shows  a  scatterplot  of  minSLP  in  the  ensemble  members 
at  1800  UTC  12  September  (i.e.,  the  initial  time  for  the 
ensemble)  versus  minSLP  at  0900  UTC  13  September 
(i.e.,  a  15-h  forecast  near  the  peak  intensity  time).  The 
strong  correlation  (—0.7)  between  the  initial  minSLP 
and  minSLP  at  the  time  of  peak  intensity  highlights  the 
importance  of  the  initial  analysis  accuracy.  The  initial 
ensemble  spread  at  1800  UTC  12  September  is  com¬ 
parable  to  or  even  smaller  than  the  typical  errors  in  the 
best-track  estimates,  which  further  demonstrates  the 
limited  predictability  of  deterministic  forecasts  even 
after  radar  observations  are  assimilated.  Thus,  the  need 
for  probabilistic/ensemble  forecasting  of  tropical  cy¬ 
clones  is  clear. 

To  further  exemplify  error  growth  between  ensemble 
members.  Fig.  14  shows  the  simulated  maximum  re¬ 
flectivity  from  one  of  the  weakest  members  and  one  of 
the  strongest  members  at  1800  UTC  12  September  and 
0900  UTC  13  September  (based  on  minSLP  at  0900 
UTC  13  September;  marked  in  Fig.  13).  Despite  having 
apparently  similar  structures  and  strengths  at  1800  UTC 
12  September,  the  two  members  diverge  tremendously 
by  0900  UTC  13  September,  again  signifying  large  un¬ 
certainties  in  the  deterministic  forecasts  of  hurricanes. 
Ongoing  studies  are  currently  investigating  both  the 
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Fig.  9.  Comparison  of  radar  reflectivity  (dBZ)  from  (left)  the  observational  composite  mosaics  (OBS),  (middle)  the  deterministic 
forecast  initialized  with  the  EnKF  analysis  at  1800  UTC  12  Sep  2007,  and  (right)  with  the  NoDA  ensemble  forecast  mean  valid  at  2100 
UTC  12  Sep,  and  0300  and  0900  UTC  13  Sep  2007,  respectively. 


mechanism  leading  to  rapid  tropical  cyclone  formation 
and  intensification  and  the  dynamics  that  lead  to  the 
rapid  error  growth  for  this  event.  Recently,  studies  re¬ 
vealed  that  upscale  growth  of  moist  convection,  such  as 
in  the  form  of  vertical  hot  towers,  may  play  a  critical 
role  in  the  internal  dynamics  (Krishnamurti  et  al.  2005; 
Montgomery  et  al.  2006).  Limited  predictability  of 
moist  convection  could  ultimately  limit  the  predict¬ 
ability  of  tropical  cyclones  (Sippel  and  Zhang  2008; 
Zhang  and  Sippel  2009),  as  is  the  case  for  extratropical 
cyclones  (Zhang  et  al.  2002,  2003,  2007;  Tan  et  al.  2004) 
or  continental  warm  season  mesoscale  convective  sys¬ 
tems  (Zhang  et  al.  2006b;  Hawblitzel  et  al.  2007;  Bei  and 
Zhang  2007). 


e.  Sensitivity  experiments 

In  the  control  EnKF  analysis  and  forecasts  discussed 
above,  we  use  the  SCL  technique  that  is  designed  to 
assimilate  dense  radar  observations  that  contain  infor¬ 
mation  about  the  state  of  the  atmosphere  at  a  wide 
range  of  scales  and  also  to  reduce  the  computational 
costs  and  sampling  errors  (refer  to  section  2e).  Four  sen¬ 
sitivity  experiments  are  performed  to  examine  the  ef¬ 
fectiveness  of  SCL  in  comparison  to  typical  covariance 
localization.  Experiments  FIXl,  FIX2,  and  FIX3  assim¬ 
ilate  the  same  SOs  as  in  the  control  EnKF  analysis  except 
that  fixed  ROIs  of  1215,  405,  and  135  km,  respectively, 
are  used  for  all  domains.  Experiment  DX30  uses  a  fixed 
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Fig.  10.  The  simulated  (a)  position,  (b)  minSLP,  and  (c) 
maxWSP  of  Humberto  from  the  1.5-km  deterministic  forecast 
initialized  with  the  EnKF  analysis  at  1800  UTC  12  Sep  2007 
(1.5KM)  in  comparison  with  the  4.5-km  forecast  and  the  best-track 
estimate. 


ROI  in  the  unit  of  horizontal  grid  spacings  in  each  do¬ 
main  (30  grid  points)  but  will  have  varying  ROIs  in 
physical  distance  (1215  km  for  Dl,  405  km  for  D2,  and 
135  km  for  D3).  DX30  is  somewhat  similar  to  the  con¬ 
trol  EnKF  analysis  with  SCL  except  that  the  ROIs  will 
not  change  in  a  given  domain. 

The  forecast  performance  of  the  WRF  simulations 
initialized  at  1800  UTC  12  September  with  analyses 
from  these  four  experiments  in  comparison  to  the  con¬ 
trol  EnKF  analysis  is  displayed  in  Fig.  15.  Neither  FIXl 
or  FIX2  with  a  fixed  ROI  captures  the  rapid  develop¬ 
ment  of  the  storm  while  DX30,  which  uses  different 
ROIs  for  different  domains,  has  decent  tropical  devel¬ 
opment,  though  it  is  significantly  weaker  than  using  the 
control  EnKF  analysis  with  SCL.  FIX3  also  simulates 
noticeably  stronger  development  of  the  storm  than  do 
either  FIXl  or  FIX2  (implying  the  benefits  and  need  for 
tighter  covariance  localization  for  the  convective-scale 
radar  observations)  but  its  results  are  progressively 
weaker  than  those  of  DX30  and  the  CNTL  analysis 
(implying  the  benefit  of  updating  larger  scales  with 
convective-scale  radar  observations  using  broader  co- 
variance  localization  in  addition  to  the  tighter  localiza¬ 
tion  for  the  convective  scales).  The  track  forecasts  from 
all  four  sensitivity  experiments  are  significantly  worse 
than  that  using  the  SCL  technique.  Nevertheless,  we 
acknowledge  that  the  assignment  of  the  numbers  of  SOs 
to  different  bins  and  the  selections  of  different  ROIs  for 
different  groups  of  SOs  are  still  empirical  at  present  and 
deserve  more  in-depth  study  in  the  future. 

Since  we  use  the  relaxation  method  as  in  Eq.  (1)  to 
inflate  the  covariance  in  order  to  avoid  filter  divergence, 
a  weighting  coefficient,  a,  set  to  0.8  (similar  to  Meng  and 
Zhang  2008a,b;  Torn  and  Hakim  2008)  implies  that  an 
overwhelmingly  large  portion  of  the  final  variance  co¬ 
mes  from  the  prior.  A  large  value  of  a  is  found  to  be 
necessary  for  imperfect-model  experiments  (Meng  and 
Zhang  2007)  to  compensate  for  unavoidable  imperfec¬ 
tions  in  the  forecast  model.  However,  in  the  presence  of 
significant  model  error  (and  with  limited  ensemble  size 
and  different  covariance  localization),  it  is  unclear 
whether  the  data  assimilation  configurations  are  opti¬ 
mal.  We  performed  two  experiments  (MIXl  and  MIX2) 
with  a  values  of  0.5  and  0.65,  respectively,  to  examine 
the  sensitivity  of  the  EnKF’s  performance  to  the  choice 
of  inflation  factor  (i.e.,  the  relaxation  coefficient  a). 
Neither  experiment  gave  better  analyses  or  forecasts 
in  terms  of  both  track  and  intensity  (Fig.  16)  compared 
to  the  CNTL  experiment  (with  weighting  coefficient 
a  =  0.8),  despite  better  agreement  between  the  mean- 
square  error  and  the  predicted  innovation  variance  (not 
shown),  which  should  be  more  desirable  (Houtekamer 
et  al.  2005).  Results  from  these  two  sensitivity  experiments 
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Fig.  11.  Comparison  of  radar  reflectivity  (dBZ)  from  (left)  observational  composite  mosaics  (OBS)  with  (middle)  those  derived  from  the 
4.5-km  control  forecast,  and  (right)  the  1.5-km  forecast  valid  at  0000  and  0600  UTC  13  Sep  2007,  respectively. 


are  consistent  with  our  recent  studies  (Meng  and  Zhang 
2008a, b)  for  regional-scale  real-data  applications.  Sev¬ 
eral  factors  may  have  contributed  to  this  mismatch  be¬ 
tween  the  mean-square  error  and  predicted  innovation 
variance:  1)  the  observational  uncertainty  is  usually  in¬ 
flated  to  be  conservative;  2)  the  consistency  assumption 
is  based  on  perfect-model  and  unbiased  statistics,  which 
are  unattainable  in  real-data  experiments;  and  3)  there 
are  some  decaying  modes  (such  as  numerical  noise)  in 
the  ensemble  variance  that  may  not  be  projected  well 
onto  the  ensemble  mean  error. 

We  also  performed  experiments  examining  the  im¬ 
pacts  of  assimilating  conventional  observations  in  com¬ 
parison  to  the  radar  observations  used  in  the  CNTL 
analysis  (Fig.  17).  Experiment  OBSl  assimilates  all  of  the 
conventional  observations  archived  in  the  Meteorologi¬ 
cal  Assimilation  Data  Ingest  System  (MADIS)  dataset 
hourly  from  0900  to  1800  UTC  12  September,  including 
radiosonde,  wind  profiler,  mesonet,  METAR  (routine 
aviation  weather  report),  maritime,  and  satellite-cloud- 
derived  wind  observations.  Forecasts  initialized  from  the 
EnKF  analysis  in  OBSl  at  1800  UTC  12  September  are 
not  capable  of  capturing  the  rapid  development  of 
Humberto  (Figs.  16c  and  16d).  On  the  other  hand,  ex¬ 


periment  OBS2  assimilates  hourly  all  conventional  ob¬ 
servations  archived  in  MADIS  along  with  the  same  radar 
observations  used  in  the  CNTL  EnKF  analysis.  Not 
surprisingly,  forecasts  initialized  with  the  EnKF  analysis 
in  OBS2  simulate  the  rapid  development  of  the  storm. 
However,  at  least  for  the  relatively  shorter  time  scales 
considered  in  this  case,  there  is  no  apparent  improvement 
in  the  analyses  and  forecasts  by  assimilating  conven¬ 
tional  observations  in  addition  to  the  radar  observations 
used  in  the  CNTL  experiment.  The  minimum  SEP  from 
OBS2  is  not  as  low  as  in  the  CNTL  run  but  it  is  well 
within  the  range  of  uncertainty  shown  in  the  control 
ensemble  forecasts  (Fig.  12).  These  additional  sensitivity 
experiments  further  demonstrate  the  importance  of  the 
convective-scale  radial  velocity  observations  provided  by 
the  coastal  WSR-88Ds.  The  success  of  the  deterministic 
forecasts  from  the  EnKF  analyses  of  CNTL  and  OBS2 
suggests  that  the  larger-scale  flow  may  be  good  enough 
not  to  adversely  impact  the  storm  development  with  the 
24-48-h  periods  that  we  examined. 

4.  Comparison  with  WRF-3DVAR 

Since  the  data  assimilation  schemes  used  in  opera¬ 
tional  forecast  models  at  NCEP  at  the  time  of  Humberto 


July  2009 


ZHANG  ET  AL. 


2119 


12Z12  18Z12  00Z13  06Z13  12Z13  18Z13  00Z14  06Z14  12Z14 


Fig.  12.  The  simulated  (a)  position,  (b)  minSLP,  and  (c) 
maxWSP  of  Humberto  by  the  ensemble  forecast  initialized  with 
the  EnKF  perturbations  at  1800  UTC  12  Sep  2007  (light  gray)  in 
comparison  to  the  deterministic  forecast  initialized  from  the  EnKF 
mean  analysis  (dark  gray)  and  the  NHC  best-track  estimate 
(black).  Analysis  and  uncertainty  until  1800  UTC  12  Sep  2007  are 
also  shown  in  dashed  curves. 


Minimum  SLP(hPa)  at  18Z12 

Fig.  13.  The  scatterplot  of  the  forecasted  minimum  SEP  at  1800 
UTC  12  Sep  (x  axis;  0  h)  and  0900  UTC  13  Sep  2007  (y  axis;  15  h) 
in  different  ensemble  members  with  the  weakest  and  strongest 
members  at  0900  UTC  13  Sep  highlighted. 


were  based  on  the  3DVAR  methodology,  here  we  use 
the  WRF  model  and  its  3D  VAR  system  to  assimilate 
exactly  the  same  observations  for  a  comparison  with  the 
EnKF  analysis  discussed  above. 

The  WRF-3DVAR  method  used  here  was  developed 
primarily  at  NCAR  and  is  now  operational  at  the  Air 
Force  Weather  Agency  (Barker  et  al.  2004).  Its  config¬ 
uration  is  based  on  an  incremental  formulation,  pro¬ 
ducing  a  multivariate  analysis  in  the  model  space.  Its 
incremental  cost  function  is  minimized  in  a  precondi¬ 
tioned  control  variable  space  where  the  errors  of  dif¬ 
ferent  control  variables  are  largely  uncorrelated.  As  in 
any  other  variational  data  assimilation  technique,  the 
structure  of  the  background  error  covariance  may  play  a 
very  important  role  in  the  3D  VAR  system.  Experiment 
“3DVAR1”  uses  the  WRF-3DVAR  default  background 
error  statistic  (its  “cv”  option  3  originated  from  an 
earlier  version  of  the  NCEP  GFS  system).  The  first 
guess  comes  from  the  9-h  WRF  4.5-km  (single,  deter¬ 
ministic)  forecast  initialized  with  the  GFS  Final  Global 
Data  Assimilation  System  (FNL)  analysis  at  0000  UTC 
12  September.  WRF/3DVAR  begins  assimilation  at 
0900  UTC  12  September  and  continues  to  assimilate 
exactly  the  same  radar  observations  as  in  the  EnKF 
analysis  hourly  until  1800  UTC  12  September,  after 
which  time  a  30-h  forecast  without  further  data  assimi¬ 
lation  is  performed.  Experiment  3DVAR2  is  performed 
in  the  same  manner  as  3DVAR1  except  that  a  30-member 
12-h  short-term  ensemble  forecast  initialized  at  0000 
UTC  12  September  (the  same  as  in  the  initial  ensemble 
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Fig.  14.  The  comparison  of  the  simulated  maximum  reflectivity  (dBZ)  derived  from  the  (left)  strongest  member  and 
(right)  weakest  member  of  the  ensemble  valid  at  1800  UTC  12  Sep  and  0900  UTC  13  Sep  2007,  respectively. 


for  the  EnKF  analysis)  is  used  to  recalculate  the  WRF/ 
3DVAR  background  error  statistics  [its  “cv”  option  5, 
as  in  Meng  and  Zhang  (2008a)]. 

Figure  18  shows  the  evolution  of  minSFP  and 
maxWSP  from  both  WRF/3DVAR  experiments  in 
comparison  to  the  corresponding  FnKF  analysis  and 
forecast.  Although  both  3D  VAR  and  FnKF  assimilate 
exactly  the  same  observations  with  the  same  model 
resolution,  both  3D  VAR  experiments  fail  to  perform 
satisfactorily.  This  result  is  in  spite  of  the  fact  that  the 
3D  VAR  analyses  may  have  a  better  fit  of  maximum 
wind  speed  at  1800  UTC  12  September  before  the  pure 
forecast  starts.  We  acknowledge  that  the  performance 
of  3D  VAR  may  be  further  improved  through  further 
tuning  of  the  background  error  covariance  (some  im¬ 
provement  is  seen  when  using  the  ensemble  to  generate 
the  background  error  covariance  in  3DVAR2)  and/or 
with  the  addition  of  initial  vortex  bogussing  (Q.  Xiao, 


NCAR,  2008,  personal  communication),  but  the  sensi¬ 
tivity  of  the  performance  of  3D  VAR  to  different  con¬ 
figurations  and  representations  of  background  error 
statistics  is  beyond  the  scope  of  this  study.  On  the  other 
hand,  there  may  still  be  room  to  improve  the  perfor¬ 
mance  of  EnKF  through  further  tuning,  which  is  also 
beyond  the  scope  of  the  current  study. 

5.  Summary  and  conclusions 

This  study  explores  the  uses  of  Doppler  radar  obser¬ 
vations  for  cloud-resolving  hurricane  analysis,  initiali¬ 
zation,  and  prediction  with  an  ensemble  Kalman  filter 
(EnKF).  The  case  studied  is  Hurricane  Humberto  (2007), 
the  first  landfalling  hurricane  in  the  United  States  since 
the  end  of  the  2005  hurricane  season  and  the  most 
rapidly  intensifying  near-landfall  storm  in  U.S.  history. 
The  storm  caused  extensive  damage  along  the  southeast 
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Fig.  15.  The  simulated  (a)  position,  (b)  minSLP,  and  (c) 
maxWSP  of  Humberto  from  four  deterministic  forecasts  initialized 
with  the  EnKF  analyses  using  fixed  radii  of  influence  (1215  km  in 
FIXl,  405  km  in  FIX2, 135  km  in  FIX3,  and  30  grid  points  for  each 
domain  in  DX30)  in  comparison  with  the  4.5-km  forecast  from  the 
control  EnKF  analysis  using  SCL  and  the  best-track  estimate. 


(a)  Track  -  ^ 


Fig.  16.  The  simulated  (a)  position,  (b)  minSLP,  and  (c) 
maxWSP  of  Humberto  from  two  deterministic  forecasts  initialized 
with  the  EnKF  analyses  using  different  covariance  relaxation  co¬ 
efficients  (a  =  0.5  for  MIXl  and  a  =  0.65  for  MIX2)  in  comparison 
to  the  4.5-km  forecast  from  the  control  EnKF  analysis  using  SCL 
and  the  best-track  estimate. 
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(a)  Minimum  SLP  (hPa)  (b)  Maximum  surface  wind  (m/s) 


Fig.  17.  The  simulated  (a)  minSLP  and  (b)  maxWSP  of  Humberto  from  the  deterministic  forecasts  initialized  with 
two  different  EnKF  analyses  (OBSl  and  OBS2)  at  1800  UTC  12  Sep  2007  in  comparison  with  the  4.5-km  forecast 
from  the  EnKF  analysis  and  the  best-track  estimate.  OBSl  assimilates  only  conventional  observations  while  OBS2 
assimilates  both  conventional  observations  and  the  radar  observations. 


Texas  coast  but  was  poorly  predicted  by  operational 
models  and  forecasters.  It  is  found  that  the  EnKF 
analysis,  after  assimilating  radial  velocity  observations 
from  three  WSR-88Ds  along  the  Gulf  coast,  closely 
represents  the  best-track  position  and  intensity  of 
Humberto.  Deterministic  forecasts  initialized  from  the 
EnKF  analysis,  despite  having  considerable  variability 
with  different  lead  times,  are  also  capable  of  predicting 
the  rapid  formation  and  intensification  of  the  hurricane. 
These  forecasts  are  superior  to  operational  forecasts, 
simulations  without  radar  data  assimilation,  and  fore¬ 
casts  initialized  with  the  assimilation  of  the  same  ob¬ 


servations  with  a  three-dimensional  variational  method 
implemented  with  the  same  forecast  model.  Moreover, 
ensemble  forecasts  initialized  with  EnKF  analysis  per¬ 
turbations  before  the  rapid  intensification  show  large 
spread  among  ensemble  members.  Such  large  spread 
further  exemplifies  the  significant  uncertainties  in  the 
deterministic  prediction  of  hurricanes,  especially  the 
intensity  forecasts. 

In  this  study,  EnKF  demonstrates  great  promise  in 
assimilating  Doppler  radar  observations  to  initialize 
hurricanes  with  detailed,  accurate  mesoscale  struc¬ 
tures.  Even  though  ground-based  radar  may  only  have 


(a)  Minimum  SLP  (hPa)  (b)  Maximum  surface  wind  (m/s) 


two  different  WRF-3DVAR  analyses  (3DVAR1  and  3DVAR2)  at  1800  UTC  12  Sep  2007  in  comparison  with  the 
4.5-km  forecast  from  the  EnKF  analysis  and  the  best-track  estimate.  The  WRF-3DVAR  default  background  error 
covariance  statistics  is  used  in  3DVAR1  while  a  12-h  ensemble  is  used  to  generate  the  background  statistics  in 
3DVAR2. 
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limited  range  in  providing  observations  for  hurricane 
prediction  with  lead  times  beyond  24  h,  this  may  be 
complemented  with  airborne  radar  with  Doppler  ob¬ 
servations  and  longer  lead  times  that  will  soon  become 
available  for  all  hurricane  aircraft  reconnaissance  mis¬ 
sions  (F.  Marks,  NOAA/HRD,  2008,  personal  com¬ 
munications).  However,  since  only  radial  winds  from 
three  radars  are  assimilated  in  this  study,  the  setup  in 
these  experiments  may  only  be  most  effective  for  short 
periods  of  time.  Beyond  a  couple  of  days,  better  con¬ 
straints  by  observations  at  the  synoptic  scales  will  be 
needed  to  further  improve  hurricane  track  and  inten¬ 
sity  forecasts. 

Future  studies  are  planned  to  examine  the  dynamics 
and  predictability  of  Humberto  with  the  EnKF  analysis 
and  forecasts.  It  remains  unclear  what  observations 
are  necessary  and  sufficient  to  define  the  initial  tropical 
cyclone  vortex  and  large-scale  environment.  Answers 
to  these  questions  have  strong  implications  related  to 
how  society  might  better  distribute  resources  to  cope 
with  future  hurricane-related  disasters.  This  is  extremely 
important  given  that  the  number  of  hurricanes  and 
their  intensity/destructiveness  are  reportedly  on  the  rise 
with  the  warming  climate  (Emanuel  2005;  Webster  et  al. 
2005). 

Despite  the  promising  performance  of  both  deter¬ 
ministic  and  probabilistic  forecasts  from  the  EnKF 
analysis,  the  intrinsic  limit  of  hurricane  predictability 
(i.e.,  in  the  face  of  nearly  perfect  observations  and  ini¬ 
tialization)  remains  unclear  in  terms  of  both  track  and 
intensity  forecasts.  The  limit  of  formation-intensity  pre¬ 
dictability  given  realistic  initial  conditions  and  model 
errors  (which  are  still  large  at  present)  in  numerical 
weather  prediction  models  may  be  alleviated  through 
improving  our  understanding  of  their  dynamics  and 
physics,  the  development  of  better  numerical  models, 
and  improved  data  coverage  and  assimilation  tech¬ 
niques.  However,  there  always  will  be  forecast  errors 
due  to  the  inherent  limit  of  predictability  arising  from 
initial  errors  with  amplitudes  far  smaller  than  any  ob¬ 
servation  or  analysis  system  (e.g.,  Zhang  and  Sippel 
2009);  these  are  errors  that  society  will  always  have  to 
cope  with  (Pielke  1997). 

Such  inherent  uncertainties  in  hurricane  forecasts 
highlight  the  need  for  developing  advanced  ensemble 
prediction  systems  to  provide  event-dependent  proba¬ 
bilistic  forecasts  and  risk  assessments.  In  practice,  de¬ 
spite  an  increasing  role  and  the  demonstrated  benefits 
of  using  ensembles  in  aiding  deterministic  hurricane 
forecasting  (Krishnamurti  et  al.  1999),  the  uncertainty 
involved  with  today’s  operational  hurricane  forecasts  is 
still  based  on  averaged  climatological  errors  and  thus  is 
not  case  dependent.  This  case  clearly  demonstrates  that 


uncertainty  can  be  quite  large  at  some  times  (e.g.,  Sippel 
and  Zhang  2008)  and  having  access  to  such  information 
in  the  operational  environment  would  serve  forecasters 
well. 
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